###########################
##################
##### Replication Rcode file for "The (still) mysterious case of agricultural protectionism"
### Authors: Quynh Nguyen, Gabriele Spilker, and Thomas Bernauer  
#### January 22, 2021


rm(list=ls())
pacman::p_load(foreign, ggplot2, plm, dplyr, reshape2, countrycode, sandwich, lmtest, MASS, 
               rworldmap, RColorBrewer, states, mice, VIM, stargazer, margins, clusterSEs, lme4, optimx,
               coefplot, survey, sampleSelection)
library(foreign)
library(cregg)
library(ggplot2)
library(prediction)
library(ggeffects)
library(dplyr)
library(tidyverse)
library(reshape)
library(ggpubr)
library(psych)
library(moderndive)
library(jtools)
library(dotwhisker)
library(stargazer)
library(cjoint)
library(gridExtra)

setwd("~/Dropbox/Agrar/Data replication") #needs to be adjusted depending on user's own workpath

## Input data 
dfagrar <- read.dta("agrardata_rep.dta") 

###################
###Prepare data for analysis
##Define variables

## Outcome variable: Agrar subsidies too low
dfagrar$farmsupport <- as.numeric(ordered(dfagrar$agrarsubtoolow,levels=c("1.Too high","2.Somewhat high",
                                                                          "3.Appropriate", "4.Somewhat low", "5.Too low")))

## Treatment: Cost frame 
dfagrar$treatment <- as.numeric(ordered(dfagrar$costframe,levels=c("0","1")))

## Covariates 
dfagrar$age <- dfagrar$Q1_recode_age
dfagrar$female <- as.numeric(ordered(dfagrar$Q2,levels=c("Male","Female")))

dfagrar$education <- NULL
dfagrar$education[dfagrar$Q4_USA_recoded=="NO COL" | dfagrar$Q4_CH_recoded=="Low"] = "Low"
dfagrar$education[dfagrar$Q4_USA_recoded=="SOME COL" | dfagrar$Q4_CH_recoded=="Medium"] = "Medium"
dfagrar$education[dfagrar$Q4_USA_recoded=="COLL+" | dfagrar$Q4_CH_recoded=="High"] = "High"
dfagrar$education <- as.numeric(ordered(dfagrar$education,levels=c("Low","Medium", "High")))

dfagrar$income <- as.numeric(ordered(dfagrar$Q22,levels=c("Less than 500 US-$",
                                                          "501 to 1,000 US-$",
                                                          "1,001 to 2,000 US-$",
                                                          "2,001 to 3,000 US-$",
                                                          "3,001 to 4,000 US-$",
                                                          "4,001 to 5,000 US-$",
                                                          "5,001 to 6,000 US-$",
                                                          "6,001 to 7,000 US-$",
                                                          "7,001 to 8,000 US-$",
                                                          "8,001 to 9,000 US-$",
                                                          "9,001 to 10,000 US-$",
                                                          "10,001 US-$ or more")))

dfagrar$employed <- as.numeric(ordered(dfagrar$Q23,levels=c("No","Yes")))
dfagrar$know_farmer <- as.numeric(ordered(dfagrar$Q26,levels=c("No","Yes")))
dfagrar$workonfarm <- as.numeric(ordered(dfagrar$Q28,levels=c("No","Yes")))

#Food security: Higher values, country should be more self-sufficient
dfagrar$prod_food <- as.numeric(ordered(dfagrar$Part3_Q5_1_scale,levels=c("Completely Disagree","Mostly Disagree", "Slightly Disagree",
                                                                          "Slightly Agree", "Mostly Agree", "Completely Agree")))  

# Traditionalism: higher values indicate higher importance of tradition
dfagrar$tradition1 <- as.numeric(ordered(dfagrar$Part3_Q6_1_scale,levels=c("Completely Disagree","Mostly Disagree", "Slightly Disagree",
                                                                           "Slightly Agree", "Mostly Agree", "Completely Agree")))                                                                           

dfagrar$tradition2 <- as.numeric(ordered(dfagrar$Part3_Q6_2_scale,levels=c("Completely Disagree","Mostly Disagree", "Slightly Disagree",
                                                                           "Slightly Agree", "Mostly Agree", "Completely Agree")))  

dfagrar$tradition3 <- as.numeric(ordered(dfagrar$Part3_Q6_3_scale,levels=c("Completely Disagree","Mostly Disagree", "Slightly Disagree",
                                                                           "Slightly Agree", "Mostly Agree", "Completely Agree")))  

dfagrar$tradition4 <- as.numeric(ordered(dfagrar$Part3_Q6_4_scale,levels=c("Completely Agree", "Mostly Agree", "Slightly Agree", 
                                                                           "Slightly Disagree", "Mostly Disagree", "Completely Disagree")))

#Generate an index of traditionalism                                                                         
dftradition <- data.frame(dfagrar$tradition1, dfagrar$tradition2, dfagrar$tradition3, dfagrar$tradition4)
alpha(dftradition) #satisfactory Cronbach alpha (.64)

dfagrar$tradition <- rowMeans(cbind(dfagrar$tradition1, dfagrar$tradition2, 
                                    dfagrar$tradition3, dfagrar$tradition4)) #using rowmeans instead of CFA


#Inequality aversion: Higher values, higher preference for equality
dfagrar$equality <- as.numeric(ordered(dfagrar$Part3_Q7,levels=c("7 We need income differences as incentive for individual effort.",
                                                                 "6", "5", "4", "3", "2", "1 Income should be distributed equally.")))

#Economic insecurity: Higher values, higher insecurity
dfagrar$econ_insecure <- as.numeric(ordered(dfagrar$Part3_Q8,levels=c("Very secure", "Secure", "Somewhat secure", 
                                                                      "Somewhat insecure", "Insecure", "Very insecure")))

#Environmental concern: Higher values, higher concern for environment
dfagrar$env_concern1 <- as.numeric(ordered(dfagrar$Part3_Q10,levels=c("Global climate change does not exist.", "… never.", "… not for many years.",
                                                                      "… in the next few years.", "... now.")))

dfagrar$env_concern2 <- as.numeric(ordered(dfagrar$Part3_Q11,levels=c("Not at all concerned", "Not too concerned", "Somewhat concerned",
                                                                      "Very concerned")))

dfagrar$env_concern3 <- as.numeric(ordered(dfagrar$Part3_Q12,levels=c("7 Economic growth and the creation of jobs should be given priority, even if this makes climate change worse.",
                                                                      "6", "5", "4", "3", "2", "1 Combating climate change should be given priority, even if this decreases economic growth and cuts jobs. ")))

#Generate an index of environmental concern                                                                       
dfenv_concern <- data.frame(dfagrar$env_concern1, dfagrar$env_concern2, dfagrar$env_concern3)
alpha(dfenv_concern) #satisfactory Cronbach alpha (.74)

dfagrar$env_concern <- rowMeans(cbind(dfagrar$env_concern1, dfagrar$env_concern2, dfagrar$env_concern3)) 


#Nationalism: Higher values, stronger nationalist sentiments 
dfagrar$nat1 <- as.numeric(ordered(dfagrar$Part3_Q14_1_scale,levels=c("Completely Disagree", "Mostly Disagree", "Slightly Disagree",
                                                                      "Slightly Agree", "Mostly Agree", "Completely Agree")))

dfagrar$nat2 <- as.numeric(ordered(dfagrar$Part3_Q14_2_scale,levels=c("Completely Disagree", "Mostly Disagree", "Slightly Disagree",
                                                                      "Slightly Agree", "Mostly Agree", "Completely Agree")))

dfagrar$nat3 <- as.numeric(ordered(dfagrar$Part3_Q14_3_scale,levels=c("Completely Agree", "Mostly Agree", "Slightly Agree", 
                                                                      "Slightly Disagree", "Mostly Disagree", "Completely Disagree"))) 

dfagrar$nat4 <- as.numeric(ordered(dfagrar$Part3_Q14_4_scale,levels=c("Completely Disagree", "Mostly Disagree", "Slightly Disagree",
                                                                      "Slightly Agree", "Mostly Agree", "Completely Agree")))

dfagrar$nat5 <- as.numeric(ordered(dfagrar$Part3_Q14_5_scale,levels=c("Completely Agree", "Mostly Agree", "Slightly Agree", 
                                                                      "Slightly Disagree", "Mostly Disagree", "Completely Disagree"))) 


#Generate an index of nationalism                                                                      
dfnat <- data.frame(dfagrar$nat1, dfagrar$nat2, dfagrar$nat3, dfagrar$nat4, dfagrar$nat5)
alpha(dfnat) #satisfactory Cronbach alpha (.7)

dfagrar$nationalism <- rowMeans(cbind(dfagrar$nat1, dfagrar$nat2, dfagrar$nat3, dfagrar$nat4, dfagrar$nat5)) 

####Distribution preference of agriculture budget 
#for environmentally friendly farming
dfagrar$environment <- dfagrar$ecofarming

#protection of small farmers 
dfagrar$farmers <- dfagrar$smallfarm

#food security 
dfagrar$food.security <- dfagrar$foodsecure

#protection of countryside 
dfagrar$country.side <- dfagrar$procountryside

#animal welfare 
dfagrar$animal.welfare <- dfagrar$animal

#food quality 
dfagrar$food.quality <- dfagrar$foodquality


######Part 2: Conjoint analysis 

#Define conjoint attributes (for US sample)

#Subset US sample which includes repondents from all conjoint tasks 
dfagrar.us <- subset(dfagrar, QCountry=="USA")

dfagrar.us$id <- dfagrar.us$respondent
dfagrar.us$choice <- dfagrar.us$accept

dfagrar.us$farmersn = rep(NA,n)
dfagrar.us$farmersn[dfagrar.us$small_farmers=="none"]="Farmers:None"
dfagrar.us$farmersn[dfagrar.us$small_farmers=="medium"]="Farmers:Medium"
dfagrar.us$farmersn[dfagrar.us$small_farmers=="high"]="Farmers:High"
dfagrar.us$Farmers <- factor(dfagrar.us$farmersn, levels=c("Farmers:None", "Farmers:Medium", "Farmers:High"))

dfagrar.us$environmentn = rep(NA,n)
dfagrar.us$environmentn[dfagrar.us$env_farming=="none"]="Environment:None"
dfagrar.us$environmentn[dfagrar.us$env_farming=="medium"]="Environment:Medium"
dfagrar.us$environmentn[dfagrar.us$env_farming=="high"]="Environment:High"
dfagrar.us$Environment <- factor(dfagrar.us$environmentn, levels=c("Environment:None", "Environment:Medium", "Environment:High"))

dfagrar.us$landscapen = rep(NA,n)
dfagrar.us$landscapen[dfagrar.us$countryside=="none"]="Landscape:None"
dfagrar.us$landscapen[dfagrar.us$countryside=="medium"]="Landscape:Medium"
dfagrar.us$landscapen[dfagrar.us$countryside=="high"]="Landscape:High"
dfagrar.us$Landscape <- factor(dfagrar.us$landscapen, levels=c("Landscape:None", "Landscape:Medium", "Landscape:High"))

dfagrar.us$foodqual = rep(NA,n)
dfagrar.us$foodqual[dfagrar.us$food_quality=="none"]="Foodquality:None"
dfagrar.us$foodqual[dfagrar.us$food_quality=="medium"]="Foodquality:Medium"
dfagrar.us$foodqual[dfagrar.us$food_quality=="high"]="Foodquality:High"
dfagrar.us$Foodquality <- factor(dfagrar.us$foodqual, levels=c("Foodquality:None", "Foodquality:Medium", "Foodquality:High"))

dfagrar.us$foodsec = rep(NA,n)
dfagrar.us$foodsec[dfagrar.us$food_security=="none"]="Foodsecurity:None"
dfagrar.us$foodsec[dfagrar.us$food_security=="medium"]="Foodsecurity:Medium"
dfagrar.us$foodsec[dfagrar.us$food_security=="high"]="Foodsecurity:High"
dfagrar.us$Foodsecurity <- factor(dfagrar.us$foodsec, levels=c("Foodsecurity:None", "Foodsecurity:Medium", "Foodsecurity:High"))

dfagrar.us$animals = rep(NA,n)
dfagrar.us$animals[dfagrar.us$animal_care=="none"]="Animals:None"
dfagrar.us$animals[dfagrar.us$animal_care=="medium"]="Animals:Medium"
dfagrar.us$animals[dfagrar.us$animal_care=="high"]="Animals:High"
dfagrar.us$Animals <- factor(dfagrar.us$animals, levels=c("Animals:None", "Animals:Medium", "Animals:High"))


#Define conjoint attributes (for Swiss sample)

#Subset CH sample which includes repondents from all conjoint tasks 
dfagrar.ch <- subset(dfagrar, QCountry=="Switzerland")

dfagrar.ch$id <- dfagrar.ch$respondent
dfagrar.ch$choice <- dfagrar.ch$accept

dfagrar.ch$farmersn = rep(NA,n)
dfagrar.ch$farmersn[dfagrar.ch$small_farmers=="none"]="Farmers:None"
dfagrar.ch$farmersn[dfagrar.ch$small_farmers=="medium"]="Farmers:Medium"
dfagrar.ch$farmersn[dfagrar.ch$small_farmers=="high"]="Farmers:High"
dfagrar.ch$Farmers <- factor(dfagrar.ch$farmersn, levels=c("Farmers:None", "Farmers:Medium", "Farmers:High"))

dfagrar.ch$environmentn = rep(NA,n)
dfagrar.ch$environmentn[dfagrar.ch$env_farming=="none"]="Environment:None"
dfagrar.ch$environmentn[dfagrar.ch$env_farming=="medium"]="Environment:Medium"
dfagrar.ch$environmentn[dfagrar.ch$env_farming=="high"]="Environment:High"
dfagrar.ch$Environment <- factor(dfagrar.ch$environmentn, levels=c("Environment:None", "Environment:Medium", "Environment:High"))

dfagrar.ch$landscapen = rep(NA,n)
dfagrar.ch$landscapen[dfagrar.ch$countryside=="none"]="Landscape:None"
dfagrar.ch$landscapen[dfagrar.ch$countryside=="medium"]="Landscape:Medium"
dfagrar.ch$landscapen[dfagrar.ch$countryside=="high"]="Landscape:High"
dfagrar.ch$Landscape <- factor(dfagrar.ch$landscapen, levels=c("Landscape:None", "Landscape:Medium", "Landscape:High"))

dfagrar.ch$foodqual = rep(NA,n)
dfagrar.ch$foodqual[dfagrar.ch$food_quality=="none"]="Foodquality:None"
dfagrar.ch$foodqual[dfagrar.ch$food_quality=="medium"]="Foodquality:Medium"
dfagrar.ch$foodqual[dfagrar.ch$food_quality=="high"]="Foodquality:High"
dfagrar.ch$Foodquality <- factor(dfagrar.ch$foodqual, levels=c("Foodquality:None", "Foodquality:Medium", "Foodquality:High"))

dfagrar.ch$foodsec = rep(NA,n)
dfagrar.ch$foodsec[dfagrar.ch$food_security=="none"]="Foodsecurity:None"
dfagrar.ch$foodsec[dfagrar.ch$food_security=="medium"]="Foodsecurity:Medium"
dfagrar.ch$foodsec[dfagrar.ch$food_security=="high"]="Foodsecurity:High"
dfagrar.ch$Foodsecurity <- factor(dfagrar.ch$foodsec, levels=c("Foodsecurity:None", "Foodsecurity:Medium", "Foodsecurity:High"))

dfagrar.ch$animals = rep(NA,n)
dfagrar.ch$animals[dfagrar.ch$animal_care=="none"]="Animals:None"
dfagrar.ch$animals[dfagrar.ch$animal_care=="medium"]="Animals:Medium"
dfagrar.ch$animals[dfagrar.ch$animal_care=="high"]="Animals:High"
dfagrar.ch$Animals <- factor(dfagrar.ch$animals, levels=c("Animals:None", "Animals:Medium", "Animals:High"))


##Figure 2: AMCEs of all conjoint attributes using choice task (US sample)

amce1 <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, id=~id)

plot(amce1)

plot(amce1) + theme(legend.title = element_blank())

plot(amce1, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank()) 


##Figure 3: AMCEs of all conjoint attributes using choice task (Swiss sample)

amce2 <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, id=~id)

plot(amce2)

plot(amce2) + theme(legend.title = element_blank())

plot(amce2, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank()) 


#######################
##### Subgroup analyses 

##Figure 4: Marginal means for all conjoint attributes by cost prime assignment (US sample)

dfagrar.us$Treatment <- factor(as.character(dfagrar.us$treatment))
levels(dfagrar.us$Treatment) <- c("Control group", "Cost treatment")

mm_cost_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                 id = ~id, estimate = "mm", by = ~Treatment)

plot(mm_cost_us, group = "Treatment", vline = 0.5)

plot(mm_cost_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure 5: Marginal means for all conjoint attributes by cost prime assignment (Swiss sample)

dfagrar.ch$Treatment <- factor(as.character(dfagrar.ch$treatment))
levels(dfagrar.ch$Treatment) <- c("Control group","Cost treatment")

mm_cost_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                 id = ~id, estimate = "mm", by = ~Treatment)

plot(mm_cost_ch, group = "Treatment", vline = 0.5)

plot(mm_cost_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold"))


######################
###Appendix 

##Figure A1: AMCEs of all conjoint attributes using rating task (US sample)

dfagrar.us$rating <- dfagrar.us$like

amce3 <- cj(dfagrar.us, rating ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, id=~id)

plot(amce3, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank()) 


##Figure A2: AMCEs of all conjoint attributes using rating task (Swiss sample)

dfagrar.ch$rating <- dfagrar.ch$like

amce4 <- cj(dfagrar.ch, rating ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, id=~id)

plot(amce4)

plot(amce4, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank()) 


##Figure A3: Marginal means for all conjoint attributes by gender (US sample)

dfagrar.us$Gender <- factor(as.character(dfagrar.us$Q2))
levels(dfagrar.us$Gender) <- c("Female","Male")

mm_gender_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                   id = ~id, estimate = "mm", by = ~Gender)

plot(mm_gender_us, group = "Gender", vline = 0.5)

plot(mm_gender_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A4: Marginal means for all conjoint attributes by gender (Swiss sample)

dfagrar.ch$Gender <- factor(as.character(dfagrar.ch$Q2))
levels(dfagrar.ch$Gender) <- c("Female","Male")

mm_gender_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                   id = ~id, estimate = "mm", by = ~Gender)

plot(mm_gender_ch, group = "Gender", vline = 0.5)

plot(mm_gender_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A5: Marginal means for all conjoint attributes by education level (US sample)

dfagrar.us$Education <- factor(as.character(dfagrar.us$education))
levels(dfagrar.us$Education) <- c("Low","Medium","High")

mm_edu_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                id = ~id, estimate = "mm", by = ~Education)

plot(mm_edu_us, group = "Education", vline = 0.5)

plot(mm_edu_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A6: Marginal means for all conjoint attributes by education level (Swiss sample)

dfagrar.ch$Education <- factor(as.character(dfagrar.ch$education))
levels(dfagrar.ch$Education) <- c("Low","Medium","High")

mm_edu_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                id = ~id, estimate = "mm", by = ~Education)

plot(mm_edu_ch, group = "Education", vline = 0.5)

plot(mm_edu_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A7: Marginal means for all conjoint attributes by income level (US sample)

##Generate 3 income levels
dfagrar.us$inccat3 <- cut(dfagrar.us$income, 3, labels=FALSE)
dfagrar.us$inccat3 <- dfagrar.us$inccat3
dfagrar.us$Income <- factor(as.character(dfagrar.us$inccat3))
levels(dfagrar.us$Income) <- c("Low","Medium", "High")

mm_inc_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                id = ~id, estimate = "mm", by = ~Income)

plot(mm_inc_us, group = "Income", vline = 0.5)

plot(mm_inc_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A8: Marginal means for all conjoint attributes by income level (Swiss sample)

dfagrar.ch$inccat3 <- cut(dfagrar.ch$income, 3, labels=FALSE)
dfagrar.ch$inccat3 <- dfagrar.ch$inccat3
dfagrar.ch$Income <- factor(as.character(dfagrar.ch$inccat3))
levels(dfagrar.ch$Income) <- c("Low","Medium", "High")

mm_inc_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                id = ~id, estimate = "mm", by = ~Income)

plot(mm_inc_ch, group = "Income", vline = 0.5)

plot(mm_inc_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A9: Marginal means for all conjoint attributes by personal tie to a farmer (US sample)

dfagrar.us$Farm_tie <- factor(as.character(dfagrar.us$know_farmer))
levels(dfagrar.us$Farm_tie) <- c("No","Yes")

mm_knowfarm_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                     id = ~id, estimate = "mm", by = ~Farm_tie)

plot(mm_knowfarm_us, group = "Farm_tie", vline = 0.5)

plot(mm_knowfarm_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A10: Marginal means for all conjoint attributes by personal tie to a farmer (Swiss sample)

dfagrar.ch$Farm_tie <- factor(as.character(dfagrar.ch$know_farmer))
levels(dfagrar.ch$Farm_tie) <- c("No","Yes")

mm_knowfarm_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                     id = ~id, estimate = "mm", by = ~Farm_tie)

plot(mm_knowfarm_ch, group = "Farm_tie", vline = 0.5)

plot(mm_knowfarm_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A11: Marginal means for all conjoint attributes by experience of working on a farm (US sample)

dfagrar.us$Farming_experience <- factor(as.character(dfagrar.us$workonfarm))
levels(dfagrar.us$Farming_experience) <- c("No","Yes")

mm_workfarm_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                     id = ~id, estimate = "mm", by = ~Farming_experience)

plot(mm_workfarm_us, group = "Farming_experience", vline = 0.5)

plot(mm_workfarm_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A12: Marginal means for all conjoint attributes by experience of working on a farm (Swiss sample)

dfagrar.ch$Farming_experience <- factor(as.character(dfagrar.ch$workonfarm))
levels(dfagrar.ch$Farming_experience) <- c("No","Yes")

mm_workfarm_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                     id = ~id, estimate = "mm", by = ~Farming_experience)

plot(mm_workfarm_ch, group = "Farming_experience", vline = 0.5)

plot(mm_workfarm_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A13: Marginal means for all conjoint attributes by traditional values (US sample)

#Generate 3 categories of traditionalism 
dfagrar.us$traditioncat3 <- cut(dfagrar.us$tradition, 3, labels=FALSE)
dfagrar.us$traditioncat3 <- dfagrar.us$traditioncat3
dfagrar.us$Traditionalism <- factor(as.character(dfagrar.us$traditioncat3))
levels(dfagrar.us$Traditionalism) <- c("Low","Mediun", "High")

mm_trad_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                 id = ~id, estimate = "mm", by = ~Traditionalism)

plot(mm_trad_us, group = "Traditionalism", vline = 0.5)

plot(mm_trad_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A14: Marginal means for all conjoint attributes by traditional values (Swiss sample)

#Generate 3 categories of traditionalism
dfagrar.ch$traditioncat3 <- cut(dfagrar.ch$tradition, 3, labels=FALSE)
dfagrar.ch$traditioncat3 <- dfagrar.ch$traditioncat3
dfagrar.ch$Traditionalism <- factor(as.character(dfagrar.ch$traditioncat3))
levels(dfagrar.ch$Traditionalism) <- c("Low","Mediun", "High")

mm_trad_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                 id = ~id, estimate = "mm", by = ~Traditionalism)

plot(mm_trad_ch, group = "Traditionalism", vline = 0.5)

plot(mm_trad_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A15: Marginal means for all conjoint attributes by inequality aversion (US sample)

##Generate 2 categories of inequality aversion
dfagrar.us$equalitycat2 <- cut(dfagrar.us$equality, 2, labels=FALSE)
dfagrar.us$equalitycat2 <- dfagrar.us$equalitycat2
dfagrar.us$Equality <- factor(as.character(dfagrar.us$equalitycat2))
levels(dfagrar.us$Equality) <- c("Low","High")

mm_equal_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                  id = ~id, estimate = "mm", by = ~Equality)

plot(mm_equal_us, group = "Equality", vline = 0.5)

plot(mm_equal_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold")) 


##Figure A16: Marginal means for all conjoint attributes by inequality aversion (Swiss sample)

#Generate 2 categories of inequality aversion
dfagrar.ch$equalitycat2 <- cut(dfagrar.ch$equality, 2, labels=FALSE)
dfagrar.ch$equalitycat2 <- dfagrar.ch$equalitycat2
dfagrar.ch$Equality <- factor(as.character(dfagrar.ch$equalitycat2))
levels(dfagrar.ch$Equality) <- c("Low","High")

mm_equal_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                  id = ~id, estimate = "mm", by = ~Equality)

plot(mm_equal_ch, group = "Equality", vline = 0.5)

plot(mm_equal_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold"))


##Figure A17: Marginal means for all conjoint attributes by economic insecurity (US sample)

#Generate 2 categories of economic insecurity
dfagrar.us$insecurecat2 <- cut(dfagrar.us$econ_insecure, 2, labels=FALSE)
dfagrar.us$insecurecat2 <- dfagrar.us$insecurecat2
dfagrar.us$Econ.Insecurity <- factor(as.character(dfagrar.us$insecurecat2))
levels(dfagrar.us$Econ.Insecurity) <- c("Low","High")

mm_insecure_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                     id = ~id, estimate = "mm", by = ~Econ.Insecurity)

plot(mm_insecure_us, group = "Econ.Insecurity", vline = 0.5)

plot(mm_insecure_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold"))


##Figure A18: Marginal means for all conjoint attributes by economic insecurity (Swiss sample)

#Generate 2 categories of economic insecurity
dfagrar.ch$insecurecat2 <- cut(dfagrar.ch$econ_insecure, 2, labels=FALSE)
dfagrar.ch$insecurecat2 <- dfagrar.ch$insecurecat2
dfagrar.ch$Econ.Insecurity <- factor(as.character(dfagrar.ch$insecurecat2))
levels(dfagrar.ch$Econ.Insecurity) <- c("Low","High")

mm_insecure_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                     id = ~id, estimate = "mm", by = ~Econ.Insecurity)

plot(mm_insecure_ch, group = "Econ.Insecurity", vline = 0.5)

plot(mm_insecure_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold"))


##Figure A19: Marginal means for all conjoint attributes by environmental concern (US sample)

#Generate 3 categories of environmental concern
dfagrar.us$envcat3 <- cut(dfagrar.us$env_concern, 3, labels=FALSE)
dfagrar.us$envcat3 <- dfagrar.us$envcat3
dfagrar.us$Env.concern <- factor(as.character(dfagrar.us$envcat3))
levels(dfagrar.us$Env.concern) <- c("Low","Medium", "High")

mm_env_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                id = ~id, estimate = "mm", by = ~Env.concern)

plot(mm_env_us, group = "Env.concern", vline = 0.5)

plot(mm_env_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold"))


##Figure A20: Marginal means for all conjoint attributes by environmental concern (Swiss sample)

#Generate 3 categories of environmental concern
dfagrar.ch$envcat3 <- cut(dfagrar.ch$env_concern, 3, labels=FALSE)
dfagrar.ch$envcat3 <- dfagrar.ch$envcat3
dfagrar.ch$Env.concern <- factor(as.character(dfagrar.ch$envcat3))
levels(dfagrar.ch$Env.concern) <- c("Low","Medium", "High")

mm_env_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                id = ~id, estimate = "mm", by = ~Env.concern)

plot(mm_env_ch, group = "Env.concern", vline = 0.5)

plot(mm_env_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold"))


##Figure A21: Marginal means for all conjoint attributes by nationalist sentiment (US sample)

##Generate 3 categories of nationalist sentiment
dfagrar.us$natcat3 <- cut(dfagrar.us$nationalism, 3, labels=FALSE)
dfagrar.us$natcat3 <- dfagrar.us$natcat3
dfagrar.us$Nationalism <- factor(as.character(dfagrar.us$natcat3))
levels(dfagrar.us$Nationalism) <- c("Low","Medium", "High")

mm_nation_us <- cj(dfagrar.us, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                   id = ~id, estimate = "mm", by = ~Nationalism)

plot(mm_nation_us, group = "Nationalism", vline = 0.5)

plot(mm_nation_us, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold"))


##Figure A22: Marginal means for all conjoint attributes by nationalist sentiment (Swiss sample)

##Generate 3 categories of nationalist sentiment
dfagrar.ch$natcat3 <- cut(dfagrar.ch$nationalism, 3, labels=FALSE)
dfagrar.ch$natcat3 <- dfagrar.ch$natcat3
dfagrar.ch$Nationalism <- factor(as.character(dfagrar.ch$natcat3))
levels(dfagrar.ch$Nationalism) <- c("Low","Medium", "High")

mm_nation_ch <- cj(dfagrar.ch, choice ~ Environment + Farmers + Landscape + Foodquality + Foodsecurity + Animals, 
                   id = ~id, estimate = "mm", by = ~Nationalism)

plot(mm_nation_ch, group = "Nationalism", vline = 0.5)

plot(mm_nation_ch, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "top") + 
  theme(strip.background = element_blank(),
        strip.text = element_text(size=8,face="bold"))


####Clear working directory and end script####
